clear all;
%% Code explanation
% 190709 use 1%Rn standard to extract Hc2, added offset for R1 and R2.

%% Define the resistance to temperature conversion equation
T_RuOx=@(R)1./exp(-1.001200439904+1.175279377873*log(R-2.2)+0.021010648114*(log(R-2.2)).^2-0.024007137972*(log(R-2.2)).^3+0.002054398372*(log(R-2.2)).^4);
%
%% Set the parameter
L_W1=64/50.6;  % sample 1: 3Sn/8PbTe
L_W2=48.5/52.9; % sample 2: 3Sn/6PbTe



%% Read the data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


root=('D:\Documents2\Rawdata\Sn\171219 MPI mK\20180223 second set of samples in plane field\'); %

filenames={'20180207_inplane_0uW_wroots';'20180208_inplane_20uW_wroots';'20180208_inplane_50uW_wroots';...
    '20180209_inplane_100uW_wroots';'20180209_inplane_150uW_wroots';'20180210_inplane_200uW_wroots'; '20180211_inplane_100uW_wrot';...
   '20180211_inplane_200uW_wrot';'20180211_inplane_400uW_wrot';'20180201_inplane_600uW_wrot';'20180213_inplane_1000uW_wrot';...
    '20180213_inplane_1200uW_wrot'};


  Rnorm=[7.6/L_W1/10 5.55/L_W2/10]'; % normal state sheet resistances of the two samples

NT=length(filenames);


colorset=jet(NT);
figure

Tmc=zeros(NT,2);
Tmc_err=zeros(NT,1);
Hc2=zeros(NT+1,2);

%% 50%Rn
ratio=0.5;
Tmax1=0.4979;
Tmax2=0.4493;

%%% 10%Rn
% ratio=0.1;
% Tmax1=0.32059;
% Tmax2=0.31182;

% %%% 1%Rn
% ratio=0.5;
% Tmax1=0.22398;
% Tmax2=0.24955;

for kk=1:NT
   
    file1=[root,filenames{kk},'_inplanescan2_waitmeas.dat'];
        fid=fopen(file1);
        data1 = textscan(fid,'%f %f %f %f %f','commentStyle','#');
        data1=data1';  
        fclose(fid);
    
     
          B=round(data1{1}*10)/10;
          R1=(data1{2}/L_W1+0.01)/10;
           R2=(data1{3}/L_W2+0.02)/10;
          T=T_RuOx(data1{4}*100);
          Tmc_err(kk,1)=3*std(T);
          theta=data1{5};
          
         B0=(min(B):0.02:max(B))';
         Bset=(min(B):0.1:max(B))';
         Bset=round(Bset*10)/10;
         RT1=zeros(length(Bset),2); % T and R1
         RT2=zeros(length(Bset),2); % T and R2
    

         for ii=1:length(Bset)
            indexB=find(B==Bset(ii)); 
            tt_err(ii)=3*std(T(indexB));
            [RT1(ii,2),j1]=min(R1(indexB));
            RT1(ii,1)=T(indexB(j1));
            [RT2(ii,2),j2]=min(R2(indexB));
            RT2(ii,1)=T(indexB(j2));
         end
          RRT1=interp1(Bset,RT1(:,2),B0,'spline');
          TT1=interp1(Bset,RT1(:,1),B0,'linear');
          

          RRT2=interp1(Bset,RT2(:,2),B0,'spline');
          TT2=interp1(Bset,RT2(:,1),B0,'linear');
          
          subplot(2,2,1),plot(Bset,RT1(:,2),'s','Color',colorset(kk,:)','MarkerSize',3);
          hold on
          subplot(2,2,3),plot(B0,RRT1,'-','Color',colorset(kk,:)','LineWidth',1);
          hold on
           subplot(2,2,2),plot(Bset,RT2(:,2),'s','Color',colorset(kk,:)','MarkerSize',3);
          hold on
          subplot(2,2,4),plot(B0,RRT2,'-','Color',colorset(kk,:)','LineWidth',1);
          hold on

          
          indexHc=find(RRT1<=ratio*Rnorm(1));
          if isempty(indexHc)~=1
              Hc2(kk,1)=B0(indexHc(length(indexHc)));
              Tmc(kk,1)=TT1(indexHc(length(indexHc)),1);
          end
          indexHc=find(RRT2<=ratio*Rnorm(2));
          if isempty(indexHc)~=1
              Hc2(kk,2)=B0(indexHc(length(indexHc)));
              Tmc(kk,2)=TT2(indexHc(length(indexHc)),1);
          end
       title={'B(T)','R(kOhm)','T(K)','averaged T(K)'};
        AA=[Bset RT2(:,2) RT2(:,1) mean(RT2(:,1))*ones(length(Bset),1)];
       
         xlswrite('FigS4B.xlsx',AA,kk);
       xlswrite('FigS4B.xlsx',title,kk);
end

subplot(2,2,1),plot(Bset,ratio*Rnorm(1)*ones(length(Bset)),'-r');
hold on
subplot(2,2,2),plot(Bset,ratio*Rnorm(2)*ones(length(Bset)),'-r');
hold on


figure
plot([Tmc(:,1)' Tmax1]',Hc2(:,1),'dr',[Tmc(:,2)' Tmax2]',Hc2(:,2),'sb');
hold on
% 
% 



